//////////////////////////////////////////////////////////////////
//                        COPYRIGHT NOTICE
// (c) Copyright 1998-2002 by Jeongnim Kim
//   All rights reserved.
//////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////
//                    CONTACT INFORMATION
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-
/*************************************************************************
 * template<class T>
 * class BoostLapackSymSolver
 *
 * Input Matrices = Lower Triangle Packed (in fortran)
 * using triangle_matrix<T,upper> 
 *
 * Expert driver to calculate M eigenstates, M < N 
 *************************************************************************/
#ifndef BOOST_NAG_SYMSOLVER_H
#define BOOST_NAG_SYMSOLVER_H

#ifndef EIGENSPACEDIMENSION_H
#include "EigSolver/EigenSpaceDimension.h"
#endif
#include <boost/numeric/ublas/config.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/traits.hpp>

// probably NAG both but just in case
#ifdef ADD_
#define f02fcf f02fcf_
#define f02fdf f02fdf_
#define f02hcf f02hcf_
#define f02hdf f02hdf_
#endif /*ADD_*/

extern "C" {
  // eigen solver Ax = lambda x A = symmetric
  void f02fcf(const char& jobz, const char& range, const char& uplo, 
	      const int& n, double* a, const int& lda,
	      const double& vl, const double& vu, 
	      const int& il, const int& iu, 
	      const int& mest, int& mconv,
	      double* w, double* z, const int& ldz, double* workd, 
	      const int& lwork, int* iwork, int& info); 
  
  // eigen solver Ax = lambda Bx A/B symmetic
  void f02fdf(const int& itype, const char& jobz, const char& uplo, 
	      const int& n, double* a, const int& lda,
	      double* b, const int& ldb, double* w, 
	      double* workd,  const int& lwork, int& ifail); 

  // eigen solver Ax = lambda x A hermitian
  void f02hcf(const char& jobz, const char& range, const char& uplo, 
	      const int& n, complex<double>* a, const int& lda,
	      const double& vl, const double& vu, 
	      const int& il, const int& iu, 
	      const int& mest, int& mconv,
	      double* w, complex<double>* z, const int& ldz, 
	      complex<double>* workd,  const int& lwork, 
	      double* rwork, int* iwork, int& info); 

  // eigen solver Ax = lambda Bx A/B symmetic
  void f02hdf(const int& itype, const char& jobz, const char& uplo, 
	      const int& n, complex<double>* a, const int& lda,
	      complex<double>* b, const int& ldb, double* w, 
	      double* rwork,
	      complex<double>* workd,  const int& lwork, int& ifail); 
}

template<class T>
class BoostNagSymSolver {

  EigenSpaceDimension& States;       

public:

  typedef T value_type;
  typedef typename boost::numeric::ublas::type_traits<T>::real_type real_type;
  typedef boost::numeric::ublas::matrix<value_type>                 HamiltonianMatrix_t;
  typedef boost::numeric::ublas::matrix<value_type>                 EigenStates_t;
  typedef boost::numeric::ublas::unbounded_array<value_type>        WorkArray_t;
  typedef boost::numeric::ublas::unbounded_array<real_type>         RealArray_t;
  typedef boost::numeric::ublas::unbounded_array<int>               IndexArray_t;
  typedef typename WorkArray_t::size_type                           size_type;

  HamiltonianMatrix_t H, S;
  EigenStates_t       Psi;  

  BoostNagSymSolver(EigenSpaceDimension& dim, bool nonorth=false): 
    States(dim), Itype(1), Range('I'), Uplo('L'), Jobz('V'){ 
    H.resize(States.N, States.N);
    if(nonortho)     
      S.resize(States.N, States.N);
  }

  ~BoostNagSymSolver() { }
  //bool getEigvectors(HamiltonianMatrix_t& H, EigenStates_t& eigvec) {
  template<class InHamiltonian_t>
  bool getEigvectors(InHamiltonian_t& h) {
    H.clear();  // clean the Hmailtonian
    copy(h,H);  // H = h*phase
    Jobz = 'V';
    Psi.resize(States.M, States.N);
    return solve(H.data().begin(), Psi.data().begin()); 
  }

  template<class InHamiltonian_t>
  bool getEigvalues(InHamiltonian_t& h) {
    H.clear();         // clean the Hmailtonian
    copy(h,H);  // H = h*phase
    Jobz = 'N';
    return solve(H.data().begin(), H.data().begin()); 
  }


  template<class InHamiltonian_t, class IterType>
  bool getEigvectors(InHamiltonian_t& h, IterType phi) {
    H.clear();         // clean the Hmailtonian
    copy(h,phi,H);  // H = h*phase
    Jobz = 'V';
    Psi.resize(States.M, States.N);
    return solve(H.data().begin(), Psi.data().begin()); 
  }

  template<class InHamiltonian_t, class IterType>
  bool getEigvalues(InHamiltonian_t& h, IterType phi) {
    H.clear();         // clean the Hmailtonian
    copy(h,phi,H);  // H = h*phase
    Jobz = 'N';
    return solve(H.data().begin(), H.data().begin()); 
  }


  template<class InHamiltonian_t>
  bool getEigvectors(InHamiltonian_t& h, InHamiltonian_t& s) {
    H.clear();         // clean the Hmailtonian
    copy(h,H);  // H = h*phase
    S.clear();         // clean the Hmailtonian
    copy(s,S);  // H = h*phase

    Jobz = 'V';
    Psi.resize(States.M, States.N);
    return solve(H.data().begin(), S.data().begin(), Psi.data().begin());
  }

  template<class InHamiltonian_t>
  bool getEigvalues(InHamiltonian_t& h, InHamiltonian_t& s) {
    H.clear();  // clean the Hmailtonian
    copy(h,H);  // H = h*phase
    S.clear();  // clean the Hmailtonian
    copy(s,S);  // H = h*phase
    Jobz = 'N';
    return solve(H.data().begin(), S.data().begin(), H.data().begin());
  }


  template<class InHamiltonian_t, class IterType>
  bool getEigvectors(InHamiltonian_t& h, InHamiltonian_t& s, IterType phi) {
    H.clear();      // clean the Hmailtonian
    copy(h,phi,H);  // H = h*phase
    S.clear();      // clean the Hmailtonian
    copy(h,phi,S);  // H = h*phase
    Jobz = 'V';
    Psi.resize(States.M, States.N);
    return solve(H.data().begin(), Psi.data().begin()); 
  }

  template<class InHamiltonian_t, class IterType>
  bool getEigvalues(InHamiltonian_t& h, InHamiltonian_t& s, IterType phi) {
    H.clear();      // clean the Hmailtonian
    copy(h,phi,H);  // H = h*phase
    S.clear();      // clean the Hmailtonian
    copy(h,phi,S);  // H = h*phase
    Jobz = 'N';
    return solve(H.data().begin(), H.data().begin()); 
  }

  inline real_type eigval(size_type i) const { return eigVal[i];}

private:

  int Itype;
  int Info;
  int Lwork, Mconv, Ifail;
  char Range, Jobz, Uplo;

  // \todo make these arrays static workspace
  IndexArray_t Iwork;
  WorkArray_t workd; 
  RealArray_t rworkd;
  RealArray_t eigVal;

  inline bool solve(double* a, double *z) {

    Mconv = 0;
    check_workspace(States.N, 64, 0);

    f02fcf(Jobz, Range, Uplo, 
	   States.N, a, States.N,
	   States.MinE, States.MaxE, States.First, States.Last, 
	   States.M, Mconv, eigVal.begin(), z, States.N,
	   workd.begin(), Lwork, Iwork.begin(), Info);

    return true;
  }

  inline bool solve(double* ap, double *bp, double* z) {

    check_workspace(States.N, 64, 0);
    f02fdf(Itype, Jobz, Uplo,
	   States.N, ap, States.N, bp, States.N, 
	   eigVal.begin(), workd.begin(), Lwork, Ifail);

    //need to copy eigen states
    if(Jobz == 'V') {
      OhmmsBlas::copy(States.N*States.M,ap,z);
    }
    return true;
  }

  inline bool solve(complex<double>* a, complex<double> *z) {

    Mconv = 0;
    check_workspace(States.N, 64, 7);
    f02hcf(Jobz, Range, Uplo, 
	   States.N, a, States.N,
	   States.MinE, States.MaxE, States.First, States.Last, 
	   States.M, Mconv, eigVal.begin(), z, States.N,
	   workd.begin(), Lwork, rworkd.begin(), Iwork.begin(), Info);
    return true;
  }

  inline bool solve(complex<double>* ap, 
		    complex<double>* bp, complex<double>* z) {

    check_workspace(States.N, 64, 3);
    f02hdf(Itype, Jobz, Uplo,
	   States.N, ap, States.N, bp, States.N, 
	   eigVal.begin(), rworkd.begin(), workd.begin(), Lwork, Ifail);

    //need to copy eigen states
    if(Jobz == 'V') {
      OhmmsBlas::copy(States.N*States.M,ap,z);
    }
    return true;
  }

  void check_workspace(int n, int workd_size, int rworkd_size) {
    Lwork = workd_size*n;
    if(workd.size() < Lwork) { workd.resize(Lwork);}
    if(Iwork.size() < 5*n) Iwork.resize(5*n);
    if(eigVal.size() < n) eigVal.resize(n);
    if(rworkd.size() < rworkd_size*n) rworkd.resize(rworkd_size*n);
  }
};

#endif

/***************************************************************************
 * $RCSfile: BoostNagSymSolver.h,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: BoostNagSymSolver.h,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/

